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ABSTRACT 

We analyze circumgalactic medium (CGM) in a suite of high-resolution cosmological re¬ 
simulations of a Milky-Way size galaxy and show that CGM properties are quite sensitive to 
details of star formation-feedback loop modelling. The simulation that produces a realistic 
late-type galaxy, fails to reproduce existing observations of the CGM. In contrast, simulation 
that does not produce a realistic galaxy has the predicted CGM in better agreement with ob¬ 
servations. This illustrates that properties of galaxies and properties of their CGM provide 
strong complementary constraints on the processes governing galaxy formation. Our simula¬ 
tions predict that column density profiles of ions are well described by an exponential function 
of projected distance d: N oc e~ d / hs . Simulations thus indicate that the sharp drop in absorber 
detections at larger distances in observations does not correspond to a “boundary” of an ion, 
but reflects the underlying steep exponential column density profile. Furthermore, we find that 
ionization energy of ions is tightly correlated with the scale height h s : h s oc £^ 4 . At z « 0, 
warm gas traced by low-ionization species (e.g., MgII and CIV) has h s ~ 0.03 — 0.07i? v ir, 
while higher ionization species (O VI and Ne VIII) have h s « 0.32 — 0.45i7 v ir- Finally, the 
scale heights of ions in our simulations evolve slower than the virial radius for 2 ^ 2, but 
similarly to the halo scale radius, r s . Thus, we suggest that the column density profiles of 
galaxies at different redshifts should be scaled by r s rather than the halo virial radius. 
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1 INTRODUCTION 


A complete picture of galaxy evolution requires an under¬ 
standing of the interplay between inflow and cooling of gas, star 
formation, and associated feedback that can drive star-forming gas 
back into the halo or beyond via large-scale winds. The current 
generation of galaxy formation models appears to capture many of 
these processes in a way that leads to galaxies with realistic prop- 


Stinson et al. 


2013 

Aumer et al. 

2013 

Hopkins et al. 2014, Agertz & Kravtsov 

2015 

Schaye et al. 

2015 

). So far, however, galaxy formation mod- 


els have been tested almost exclusively against properties of the 
stellar component of galaxies (with some recent exceptions, e.g., 
|Hummels et al.|2013[|Fbrd et al.|2013||20l5||Schaye et al.|2015) . 

Yet, we know that stars comprise only a relatively small frac¬ 
tion of the baryon budget. A variety of observational constraints 
and inferences from the abundance matching technique (e.g., see 
Fig. 10 and 11 in [Kravtsov et al.||2014[ for a compilation of re¬ 
cent constraints) shows that even the galaxies that are most effi¬ 
cient in forming stars (forming in halos of Mh ~ 10 12 M©) con¬ 
vert at most « 30 — 40% of the available baryon budget into stars. 
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The bulk of the remainder of the baryon budget is outside of the 
disk, because cold gas fractions in the evolved galaxies we ob¬ 
serve today are also small. The “missing” baryons are thought to 
reside in the the circumgalactic medium (CGM) around galaxies, 
which can be thought of as a halo of tenuous gas resulting from the 
complex interactions between static gaseous halo, inflowing gas, 
and feedback-driven outflows from galaxies. Studying the galactic 
baryons in the CGM is therefore paramount for the understanding 
of the feedback processes that regulate inflows and outflows of gas. 


Although details of such processes are still debated, efficient 
feedback recipes have been implemented in galaxy formation sim- 


ulations (e.g., Hopkins et al. 

2011, Brook et al. 2012, Stinson 

et al. 2013, Yogelsberger et a 

2013; Agertz et al. 2013, Roskar 


et al.||2014[ |Trujillo-Gomez et aI7||201 5| >. At the low mass end 

(Mh < 10 12 M 0 ), stellar feedback in the form of energy and mo¬ 
mentum injection from radiation pressure, stellar winds, and super¬ 
nova explosions is believed to be responsible for the suppression of 
star formation ( [Dekel & Silk|1986l|Efstathiou|2000| ). On the other 
hand, active galactic nuclei (AGN) feedback is thought to be crit¬ 
ical in limiting star formation and stellar masses of galaxie s at the 


high mass end (Mh > 10 M 0 ; e.g., Silk & Rees 1998 


Benson 


et al. 2003). More recently, cosmic ray (CR) driven winds have also 
been shown to be a promising feedback mechanism in the context 
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of galaxy formation simulations (e.g., Socrates et al. 2008, Booth] 
|et al.|2013l|Salem et al.|2014) . 

It is quite likely that stellar feedback leaves a specific imprint 
on the properties of the circumgalactic gaseous halos (e.g., |Barai] 
|et al.|2013l|Hummels et al.|2013||Suresh et al.|2015| >. In particu¬ 
lar, comparisons between simulated CGM and observations have 
indicated that the observed incidence of absorbers in the CGM fa¬ 
vors larger intensity of feedback processes (e.g., |Stinson et al.|2012[ 
|Hummels et al.|2013l|Suresh et al.|2015) . 

Observations of the CGM have advanced tremendously in re¬ 
cent years. The Cosmic Origins Spectrograph (COS; [Green et al.| 
|2012) on board the Hubble Space Telescope (HST) has offered a 
unique opportunity to systematically probe low-redshift (z ~ 0 ) 
galactic halos in the UV regime. Commonly observed absorp¬ 
tion lines arising in the CGM are the Lya A1215, Sill A1260, 
Si III A1206, CII A1334, CIV AA1548, 1550 (e.g., |Chen|[20l2l 
|Borthakur et al.|20T3l|Liang & Chen|2014[|Bordoloi et al.|2014|), 

and OYI AA1031, 1036 (e.g., |Prochaska et al.||2011| |Tumlinson| 

|et al.| 201 1; Johnson et al. 2015J. In the intermediate redshift range, 

many researchers have probed the gaseous halos using optical 
ground-based facilities via Mg II ( |Chen et ak| 2010, Gauthier et al.| 
|2010|[Bordoloi et al.|2011| ). Similar studies for the UV transitions 
have also been conducted in the high redshift universe (z & 2 ; 
|Steidel et al.|2010||Rudie et al.|2012|[Turner et al.|2014| ). 

In this paper, we explore how the choices of parameters con¬ 
trolling star formation and stellar feedback in cosmological galaxy 
formation simulations affect the properties of the CGM around sim¬ 
ulated galaxies, such as the profile of the absorber column densities. 
We also use the simulation in which the CGM properties are closest 
to observations to explore the physical origin of the CGM absorbers 
and the processes that shape their properties. 

The paper is structured as follows. In section 2, we outline the 
simulation suite used in our study, including details of star forma¬ 
tion and the stellar and cosmic ray feedback models. In section 3, 
we describe the post-processing of simulations and our CGM anal¬ 
yses. Main qualitative features of the CGM in our simulations are 
discussed in section 4, where we demonstrate that the radial profiles 
of the CGM absorber column density have an exponential form and 
evolve only very weakly with redshift. In addition, we show that 
profiles of the CGM spanning four decades of stellar mass and 11 
billion years in cosmic time can be cast as a single profile using a 
simple re-scaling of radius (d d/r s ). We then present compar¬ 

isons between our predictions from simulations and observations of 
the CGM in section 5, where we show that observations also appear 
to obey scaling with halo scale radius. In section 6 we discuss our 
results and their implications, as well as compare them to results of 
other recent studies of the CGM in simulations. Finally, in section 
7, we summarize our results and conclusions. 


2 GALAXY FORMATION SIMULATIONS 

In this study, we use a suite of galaxy formation simulations, 
which consists of “zoom-in” re-simulations of the evolution of a 

Milky Way-sized halo started from the same initial conditions, but 
with different parameters of star formation and stellar feedback 
recipes, as described in Age rtz & Kravtsov] ( |2015| >. In addition, 
we use a new re-simulation with a cosmic ray feedback model de¬ 
scribed in |Booth et al.[ ( |2013] ). 

All simulations are run using the Adaptive Mesh Refinement 
(AMR) code RAMSES ( |Teyssier|2002| ) from the same initial condi¬ 
tions in the WMAP5 ACDM cosmology with Qa = 0.73, n m = 


0.27, Q b = 0.045, cr 8 = 0.8 and H 0 = 70kms“ 1 Mpc“ 1 |k^| 
jmatsu et al.|2009| ). At z = 0, the halo masses of the simulated ~ L* 
galaxies have similar values of M v i r ~ 10 12 M©, corresponding to 
a virial radius of P v ir ~ 260 kpc. The dark matter particle mass 
in the highest resolution region is ttidm = 3.2 x 10 5 M Q . The 
physical resolution at maximum refinement level reaches Ax ~ 75 
pc. A summary of the parameters of our simulations is provided in 
Table 1. 

Unless otherwise stated, all quantities in the paper are given in 
physical units. We briefly outline the main ingredients of the sim¬ 
ulations in the following subsections. A more detailed description 
can be found in Agertz & Kravtsov (2015), while the specifics on 
the implementation of stellar and cosmic ray feedback are detailed 
in | Agertz et al.| ( |2013] ) and |Booth et al.| ( [20T3] ), respectively. 

2.1 Star Formation 

In each star forming cell, the number of star particles N is 
determined from the Poisson distribution P(N\Xp) with the mean 
A p — p*Ax 3 At/m*, where star formation rate is given by: 

P* = fa 2 ^es- ( 1 ) 

iff 

Here /h 2 is the local mass fraction of molecular hydrogen com¬ 
puted using the model of Krumholz et al. ( 2009), as described in 
Agertz & Kr avtsov|(|2015|, p g is the local gas density of a given 
cell, tf£ — y / 37r/(32Gp g ) is the free-fall time, and is the star 
formation efficiency per free-fall time assumed to be constant in 
time and space. The simulations analyzed in this paper use two 
constant values of the efficiency: es — 1% and 10%. The param¬ 
eter ra*, set to 10 4 M© in our simulations, is a unit mass of star 
particles and all star particles are created as a multiple of it ( |Rasera| 
|& Teyssier|2006[[Dubois & Teyssier|2008| >. 

Parameters of the star formation recipe have a direct impact 
on the CGM because the feedback strength depends on the mass 
and spatial distribution of young stars, as we will discuss in section 
[4](see also | Agertz & Kravtsov|2015[ [Stinson et al.|2012| . 

2.2 Stellar Feedback 

Our simulations account for energy, momentum, mass, and the 
injection of heavy elements (with an effective yield of 1% - 3%) 
from type la supernovae (SNIa), type II supernovae (SNII), stellar 
winds and radiation pressure from massive stars, as well as for the 
secular mass loss into the surrounding interstellar medium (ISM) 
from the AGB stars. Injection of momentum, energy and mass is 
done continuously during each simulation time step, which is com¬ 
puted according to the Courant-Friedrichs-Lewy (CFL) condition. 
In the regions of active star formation and feedback the time steps 
become as short as a few thousand years. 

When a supernova explodes, we retain a fraction, /fb, of 
the explosion energy in a separate energy variable En>, which is 
evolved using an equation similar to that of the internal energy of 
the gas, but with the explicitly specified dissipation time scale, fdis* 

4(Sfb) + V-(Sfb«gas) = -PfbV-Vgas- (2) 

Cft tdis 

This energy variable is used to define an additional contribu¬ 
tion to gas pressure, Pfb = (7 — 1 )£’«,, which can be considered as 
a crude approximation for effective pressure of the hot gas < [Agertz| 
|& Kr avtsov 2015) or subgrid turbulence unresolved in simulations 
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Table 1 . Simulation suite 


Simulations 


Description 

Final Redshift 

KMT09 models, feedback energy variable Er,, /r> 

= 0.5, tdis = lOMyr 



ALL_Efb_e010 

ALL_Efb_e001 

ALL_Efb_e001_5ESN 


All feedback processes, eff = 10% 

All feedback processes, eff = 1% 

All feedback processes, Esnii = 5 x 10 51 erg, eff = 1% 

z = 0 

z = 1.5 

z = 0 

Cosmic rays model, feedback energy variable Ecu 

= £crT?snii 



ALL_e010_CR 


All feedback processes, £cr = 10%, eff = 10% 

is* 

II 

b 


and can allow for more efficient transfer of supernova momentum 
and energy into the surrounding ISM. The remaining 1 — /fb frac¬ 
tion of energy is added to the thermal energy of the gas in the 
cell containing the stellar particles. The dissipation time scale is 
adopted to be tdis = 10 Myr for all simulations, comparable to a 
couple crossing times in massive giant molecular clouds. 

In addition, we model momentum injection due to radiation 
pressure, stellar winds, and supernovae by directly adding momen¬ 
tum into surrounding cells using the rate: 


Prad = (m + mrm) 



(3) 


where tir is the optical depth in the infrared (IR) band, L(t) 
is the luminosity of the stellar population, taken from the stellar 
evolution code STARBURST99 ( [Leitherer et al.|jT999| ). The first 
term corresponds to direct momentum injection from UV photons, 
where rji is chosen to be unity. The second term corresponds to 
momentum transfer by IR photons re-emitted by dust after the 
UV photons are absorbed, where 772 = 2 is adopted in |Agertz &] 
|Kravtsov] ( |2015| ). 


2.3 Cosmic Ray Feedback 

Simulations of isolated disk galaxies in |Booth et al.| ( [2013} and 
Salem et al. ( 2014) show that cosmic ray (CR) driven winds contain 
more cool/warm gas ( T < 10 5 K) than the winds driven via direct 
momentum and thermal energy injection from supernovae. This is 
because in the CR-driven winds the gas is gradually accelerated 
by a large-scale pressure gradient established by the cosmic rays 
diffusing out of the disk, while in the standard winds produced by 
energy and momentum injection at SN sites the gas is launched at 
large velocity and is then shock-heated to high temperatures. There¬ 
fore, we also consider a new re-simulation of the same halo as in the 
other simulations, but instead incorporates the CR feedback using 
the isotropic diffusion model described in |Booth et al.| ( [2013) . 

Briefly, the CRs are assumed to be produced by supernovae 
with a fraction £cr of the SN energy converted into the energy of 
cosmic rays: Ecu = £cr^sn- The rest of Eq X , is added to the 
thermal energy of the gas in the SN parent cell. In the simulation 
used in this analysis, we adopted £cr = 10 % consistent both with 
recent empirical evidence for CR acceleration in the SN remnants 
(|Ackermann et al.|2013|> and theore tical models (Maik ov" & Drury| 
|2001||Caprioli & Spitkovsky|2014| ). 

The cosmic rays are modelled as an ultra-relativistic ideal fluid 
with 7 cr = 4/3, which exerts pressure Pen = (tcr — 1)E?cr, 
where Ecu is the energy density of the cosmic ray particles fol¬ 


lowed using a separate variable. The variable is evolved using the 
advection-diffusion equation, which includes terms corresponding 
to cooling losses due to decays and Coulomb interactions of CR 
with gas and magnetic fields, an isotropic diffusion term assuming 
the constant diffusion coefficient offt = 3xl0 -27 cm 2 s -1 , and a 
source term due to the SN remnant production of cosmic rays. We 
also include the heating term into the energy equation of baryon 
gas due to interactions with CRs (see details in |B ooth et al. 12013] ). 

It is important to keep in mind that properties of galaxy form¬ 
ing in a halo and properties of its CGM may depend on the choices 
of our model parameters (e.g, the diffusion coefficient). The param¬ 
eters have sizeable uncertainties and thus, in principle, they need to 
be varied to explore their full effect on the CGM (e.g [Salem et al.| 
|2014| >. In this pilot study, we choose to investigate results only for 
the fiducial parameter values, chosen to be consistent with current 
observational constraints and theoretical expectations. 

Note also that CRs can lead to gamma ray emission via pion 
production. However, in our model, we find that only < 1% of 
the total energy would be emitted as gamma rays, consistent with 
a detailed CR modelling for the Milky Way ( [Strong et al.||2010| ). 
Therefore, the gamma rays should not have a significant effect on 
the low-density circumgalactic medium, but can be used as inter¬ 
esting observational constraints on the cosmic ray feedback in the 
future. 


3 CGM ANALYSIS 

3.1 Lines of sight and ion column densities 

Observed quasar spectra reveal existence of a wide range of 
ion species tracing different temperatures, which indicates that the 
CGM is multi-phase. It is thus important to use the entire range 
of ion species to test whether simulations reproduce the observed 
CGM structure. 

We compute abundances of all the ions for which extensive 
observational measurements have been reported in the literature. 
These ions cover a wide range of ionization energies, Ei on , as 
shown in Table 2. We will adopt a convention to name species with 
E ion < 54.4 eV as low ions (e.g., HI, Mg II, Sill, Si III, Si IV 
and CII), those with 54.4 < .Eion < 100 eV as intermediate ions 
(CIV), and those with Ebon > 100 eV as high ions (O VI). 

Following |Smith et ak] ( |2011| ) and |Hummels et al.| ( |2013| >, we 
assume that the number density of ion j of an element X can be 
computed as: 

n Xj (n H ,T, Z) = f Xj (T,n H )fx(Z)n H , (4) 

where f x (Z) — n x /riH is the fraction of the element X relative 
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z 


Figure 1. Star formation histories of the main galaxy in our simula¬ 
tions compared to the average histories f or observed galaxi es in halos of 
M v i r (z = 0) = 10 12 M o derived by Behroozi et al~| |2013| with dark and 
light gray shaded bands showing the one- two-sigma confidence regions of 
their constraints, respectively). The SFHs in simulations are averaged over 
time bins of AtgF = 100 Myr. 


to hydrogen and fx 6 (T, tih) = riXj/nx is the ionization fraction 
of X in state j. 

We compute an interpolation table of fx i (T, uh, Juv) using 
the CLOUDY code (version 13.02; last reviewed in |Ferland et al.| 
|2013] > under assumption of ionization equilibrium for a grid of 
riH, T, and redshift (z, corresponding to different cosmic mean 
Juv (z)). For Juv we adopt the redshift dependent ionization back¬ 
ground HM05 (Ha ardt & Ma dau 2012), which includes contri¬ 
bution from QSOs and galaxies and assume that gas is optically 
thin throughout the CGM. We assume solar pattern of heavy el¬ 
ement abundances ( |Asplund et al.||2009) and solar metallicity of 
Zq = 0.02, and approximate the mean molecular weight by a con¬ 
stant /I — 0.62 for ionized gas with T > 10 4 K, and fi — 1 
otherwise. 

The column density of ion Xj along a given line of sight 
through simulation volume is computed as: 

Nxj = J nXj(riH,T, Z, z)dl, (5) 

where rtXj is 3D number density. 

To sample the CGM uniformly under all possible viewing an¬ 
gles, we uniformly sample impact parameter d between the putative 
QSO and simulated galaxy and solid angle at the point in the line 
of sight (LOS) closest to the galaxy. To make sure that column den¬ 
sities for the LOSs near the edge of the box (i.e., large d) are not 
biased low, we have chosen a large box size Lbox = 1 Mpc (co¬ 
moving) compared to the virial radius of the simulated dark matter 
halo. Furthermore, to ensure that column densities are not biased 
up to d ~ 2.5i^vir, we do not include the LOS with path lengths 
0.27Lbox in our analysis and plots. 


All feedback + Efb,eff = 10% All feedback + Efb,eff = 1%,5 x Esnii 



Figure 2. Maps of the total gas column density, metallicity, and mass- 
weighted temperature of the fiducial and 5 -Esnii models at z = 0. The 
maps show the outflow streams aided by the boosted supernova events while 
the fiducial run maintains a normal star-forming disk. As a consequence, gas 
metallicity is well-mixed in 5L?snii model out to large radii at 10% solar 
to solar value. In contrast, despite the disk in the fiducial run having so¬ 
lar metallicity, the gas hovers at a few percent solar in the circumgalactic 
space. Note that the color of the star-forming disk in the fiducial run has 
been clipped in order to ensure the dynamic range of the color map in the 
halo. 


3.2 Curve of Growth Analysis 

In observational samples, direct column density measure¬ 
ments are often available only in relatively high signal-to-noise and 
high resolution spectra, while some studies only provide the equiv¬ 
alent width W r of absorbers. To combine all the data sets in a uni¬ 
form manner, we follow [Hummels et al. 1(2013) and apply a curve 
of growth analysis to convert from column density to equivalent 
width and vice versa using an approximation of equation (9.8) & 
(9.27) in [Praine|pOTl| : 



where 
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W r is the rest frame equivalent width in units of wavelength, to is 
the optical depth given the column density Ni and b parameter of 
the specie in state Z, 7 i u , f\ u and Ai u are the intrinsic width, oscilla¬ 
tor strength and wavelength of the specie transition from state l to 
u. The atomic data for all transitions are taken from a compilation 
study by |Morton] ( |2003| ). 

To compute N from W r , we invert the function W r (N, b) 
using its bivariate spline approximation. When only W r measure¬ 
ments are available, the set of b values are determined by a Gaus¬ 
sian width <r, as b = \/2or, which are taken from the literature 
(e.g., see |Hummels et al.|2013||Liang & Chen|2014| >. To estimate 
W r from N in simulations, we treat b as a free parameter, because 
simulations do not model non-thermal components of the gas. We 
discuss this in detail in section[5] 


4 CGM IN SIMULATIONS 

In this section we discuss the main qualitative features of the 
CGM in our simulated galaxies. We discuss their dependence on 
the feedback models (different implemented physics), as well as 
varying parameters of star formation and feedback recipes. 

Figure [T] shows star formation histories of all the simulations 
used in this study compared to the range of star formation his¬ 
tories derived for observed galaxies hosted by halos of similar 
mass Behroozi et al.|p0l3) . We can see that the main progeni¬ 
tor in the ALL_Efb_eO 01 run has star formation rates at z > 2 
that are up to an order of magnitude higher than the rate inferred 
by |Behroozi et al.| ( |2013| ). This leads to an overestimate of stel¬ 
lar mass of the final object by a factor of two and creation of a 
dense central spheroid in the stellar distribution of the galaxy at 
low z. In other runs, feedback is sufficiently efficient to suppress 
star formation to the level consistent with the inference of B ehrooz i 
|et al.H2013| . At lower z the star formation rate in the fiducial run 
ALL_Efb_eO 10 continues to track the empirical SFH, while the 
rates in the ALL_Efb_eO01_5ESN and ALL_e010_CR runs fall 
below the gray band. As we discuss below, the bursty nature of 
star formation in the ALL_Efb_e001_5ESN simulation, however, 
allows it to drive significant outflows out of the central disk even 
at low redshifts, as can be seen in Figure [2] The ALL_e010_CR 
run tracks the SFH derived from observations quite well down to 
z — 1.5, although the SFH starts to decrease sharply after this 
epoch. 

As shown by |Agertz & Kravtsov] < |2015[ |2016| ), the run 
ALL_Efb_eO 10 not only reproduces the observed star formation 
history (SFH), but produces a galaxy with realistic stellar mass- 
halo mass ratio, bulge-to-disk ratio, metallicity, and the Kennicutt- 
Schmidt relation. On the other hand, the ALL_Efb_eO01_5ESN 
run, with less efficient star formation but more energetic SNe 
events, under-produces stellar mass at low z, and fails to produce a 
disk-dominated galaxy. 

4.1 Main features of the CGM in simulations 

Figures [2] [5] and [4] show the overall thermodynamic prop¬ 
erties of the CGM and column density maps of specific ions 
at z — 0. Here we focus on the runs ALL_Efb_e010 and 
ALL_Efb_eO01_5ESN and will compare results with other runs 
in the next subsection. The figures clearly show that these two 
runs produce not only very different central galaxies, but very dif¬ 
ferent CGMs. The CGM of the ALL_Efb_eO 10 run is hot, low- 
metallicity, and is almost devoid of gas with T ~ few x 10 5 K in its 


inner regions. In contrast, the CGM in the ALL_Efb_eO01_5ESN 
simulation does have such warm gas due to gas lifted out of the 
disk by recent outflows. These can be seen as “tails” and “streaks” 
in the column density map. 

Figures [3] and [4] also show that the radial distribution of the 
high ions is more extended compared to the low ions. Overall, anal¬ 
yses of |Ford et al.| ( |2014] ) and |Ford et al.| ( |20l3] ) show that higher 
ionization energy ions like, O VI, originate from the gas that was 
ejected from the disk with outflows at higher redshifts, while low- 
ionization energy ions originate in gas associated with recent out¬ 
flows. 

4.2 Cosmic Ray Feedback on the CGM 

Recent studies (e.g., |Booth et al.| [2013] |Salem & Bryan| 
|2014| ) have shown that CR-driven winds can create large out¬ 
flow mass loading factors, especially in dwarf galaxies. Fig¬ 
ures \5\ and [6] show that both boosted SN-driven feedback in the 
run ALL_Efb_e001_5ESN and CR-driven outflows in the run 
ALLeOlCLCR produce column density profiles of similar shape 
and extent. 

However, the processes that launch gas out into the CGM in 
these simulations are quite different. In the ALL_E f b_e 0 01 _5E SN 
run outflows are driven by energy and momentum injection at the 
sites of the SNe explosions in the disk and are quite bursty. In the 
CR run, on the other hand, the outflows are driven by a large-scale 
pressure gradient established by the cosmic rays diffusing out of 
the disk. The acceleration to this gradient is more gradual, which 
results in lower outflow velocities and resulting lower gas temper¬ 
ature (see |Booth et aL]|2013| |Salem & Bryan||20T4] | Salem et~ah| 
2014 for more details on the wind properties in the CR feedback 
model). The resulting distribution of the CGM in the CR run is 
much smoother than in the ALL_Efb_e0 01_5ESN run, as can be 
seen in the right column of Figure [5] 

4.3 Evolution of the CGM in runs with different feedback 
models 

Figure [5] also shows evolution of the CGM in our different 
runs at different z. We choose to show Mg II AA2796, 2803 lines 
as an illustrative example, but similar behavior is seen for other 
ions that trace the cold/warm gas (e.g., CIV). Figure [5] shows that 
the extended gaseous halo is established by z — 3 in all runs ex¬ 
cept ALL_Efb_e001. In the case of ALL_Efb_e001, the gas in 
the halo collapses onto the central disk where it is converted into 
stars, which results in high star formation rates (see Fig. [TJ. The 
ubiquity of the extended CGM in all runs at high z is due to the 
shallow gravitational potential at early epochs, which allows even 
moderate feedback to eject gas from the disk and launch it far into 
the halo (and even beyond). It is thus more difficult to use CGM for 
differentiating feedback recipes at high redshift (z ^ 2 — 3). 

At redshifts z ^ 2, the galaxy in the simulation 

ALL_Efb_e010 continues to form stars with the rate that matches 
the empirical average SFH of Beh roozi et al.| ( |2013| ) for galaxies 
of this mass. However, the feedback implementation fails to drive 
winds far into the halo at z < 1; the SN action at these red¬ 
shifts is limited to small fountains and stirring of the gas within 
the disk. On the other hand, the runs ALL_Efb_e001_5ESN and 
ALLe_010_CR drive significant outflows down to the smallest red¬ 
shifts to which they were run. The winds lift cool gas from the disk 
into the halo in extended plumes, which are prominent in Figures 
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Figure 3. The spatial column density distribution (500 kpc across; R v i r ~ 250 kpc) of the chemically enriched CGM as probed by HI and Mg II at z = 0.0. 
Solid circle represents R v i r , dashed circle is 4 r s and dotted circle is r s . Our fiducial model with star-formation efficiency eff = 10% produces a stellar 
component consistent with observations (Agertz & Kravtsov|2015|[20T6) but fails to produce an extended multi-phase gaseous halo. On the other hand, our 
strong feedback model (with five times supernova energy, 5 x Esnii) matches the CGM profiles despite destroying the star-forming disk. This shows that the 
CGM provides sensitive orthogonal constraints on galaxy formation, especially feedback recipes. 


Another interesting feature of the evolution shown in Figure[5] 
is that the extent of the high-column density area increases visibly 
slower than the virial radius (shown by the solid lines). The same 
can be seen in the column density profiles at different z in Figure 
[6] in which profiles are rescaled by host halo R v i r (z) at the cor¬ 
responding redshift. For example, although visible physical extent 
of the Mg II distribution increases from z = 3 to z = 1 in runs 
ALL_Efb_e001_5ESN and ALL_e010_CR in Figure [5] the extent 
in units of virial radius in Figure[6]decreases over the same redshift 
interval. 


The virial radius of host halos has now been used in a number 
of studies to rescale the input parameter of galaxies at different red- 
shifts in order to compare their pr ofiles ([C hurchill e t al.|2013| [ Werk 

|et al.|2014] |Liang & Chen|2014[ |Bordoloi et al.|2014) . However, 

the evolution of the profile in simulations indicates that the actual 
scaling may be different. This motivates a closer look at the evolu¬ 
tion of the column density distribution and corresponding scaling 
across redshifts. 


4.4 Evolution of the radial distribution of the CGM column 
density and rescaling with halo scale radius 

The boundary of dark matter halos is usually defined as the 
radius enclosing a given density contrast A (e.g., A= 200) rela¬ 
tive to a reference density, such as the critical, p C rit(^), or mean, 
Pm(z), density of the universe. Common choices are R v i r (adopted 
in this paper), i? 2 oom and i? 200 c- These definitions of the radii are 
loosely motivated by simple models of spherical top-hat collapse. 
For real halos the increase of the virial radius roughly tracks the 
actual physical splashback radius R sp ( [More et al.||2015| , which 
delineates the outer region enclosing recently accreted matter. Dur¬ 
ing epochs when halos accrete quickly, the entire halo profile scales 
well with the virial and splashback radii. However, when the accre¬ 
tion rate slows down the inner density profile remains almost static 
in physical units (e.g., |Cuesta et al.|2008[ ) and thus no longer scales 
well with the outer virial or splashback radius. Galaxies at z ^ 2 
are in the latter regime, which means that scaling with the virial 
radius may not be optimal, especially since observations probe the 
CGM mostly in the inner halo (^ 4 — 6r s ). 
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Figure 4. Similar to Figure[3] but for intermediate ion CIV and high ion O VI at z = 0.0. The multi-phase gaseous halos are progressively more extended 
from low and intermediate ions (e.g., Mg II and CIV) to high ion (O VI). 


Evolution of the CGM for run ALL_Efb_e001_5ESN is il¬ 
lustrated in Figure [7] which shows evolution of the column den¬ 
sity profiles of different ions from z = 3 to z = 0 in physical 
kpc and scaled with ,R v i r in the left and middle columns. Over the 
range of redshifts shown the dark matter halo mass grows from 
M v ir = 3 x 10 11 M© to M v ir = 10 12 M©, while R v i r grows by 
by a factor of 4.6. The figure shows that the profile in physical kpc 
evolves only very weakly over 11 billion years since z = 3 and that 
scaling with virial radius introduces spurious evolution. 

Recently, [More et ak] < |2015) advocated a mass and radius defi¬ 
nition to characterize evolution of the inner regions of halo profiles 
based on halo scale radius r s , defined as the radius where logarith¬ 
mic slope of the dark matter density profile is —2. In particular, 
they argued that the virial radius in the fast accretion regime is ap¬ 
proximately equal to 4 r s , while in the slow accretion regime the r s 
and 4 r s cease evolution ( [Bullock et al.|200f| ) while R v i r continues 
to grow due to pseudo-evolution ( [Diemer et al.||2013| ), as can be 
seen in Figure [8] 

If the gravitational potential well in the inner region plays a 
major role in controlling the extent of the CGM, it is natural to 
assume that the CGM profiles should scale with r s . Indeed, the 


right column of Figure [7] shows that the column density profiles 
scale with r s much better than with i^virQ 

Note that the scaling of profiles around galaxies of different 
stellar masses at a given z with r s (or any multiple of it) is very 
similar to scaling with R v i r at the same z because r s — i2 v ir/c v ir 
and c v ir oc M^ 0-1 is a weak function of mass. However, both R v i r 
and c v ir evolve fast with redshift in a way that leaves r s almost 
constant at z < 2 for galaxy sized halos. Thus, re-scaling of profiles 
of objects at different z with r s and R v i r is quite different. 

Although concentration at a fixed halo mass exhibits scatter 
of « 0.14 dex, which will be added in quadrature to the scatter in 
the M* at fixed M vir during conversion from M* to M v i r , Figure[7] 
shows that such additional scatter may be worth the reduced bias in 
the redshift rescaling. As we will show in the next section, observed 
column density profiles also appear to favor rescaling with r s . 


1 To estimate the scale radius r s = R v i r /c v i r , we use the 
halo-concentration model of Die mer & Kravtsov| {2015), im¬ 
plemented in the publicly available python code Colossus: 
http://bdiemer.bitbucket.org 
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All feedback + E^, eg = 10% All feedback + Efb, eg = 1% All feedback + Efb, eg = 1%, 5 x Esnii All feedback + Efb,£cR — 10% 



Figure 5. Model predictions of Mg II column density maps for different redshifts. The region is 250 kpc (physical) on a side. Solid circles represents R v - 1T (it 
is outside the region at z = 0) and dashed circle represents 4 r s . At high redshifts, 4 r s tracks the growth of R v i r , but at z < 2 R v [ r becomes larger and grows 
faster than 4 r s . For models that match observations of the CGM (ALL_Efb_e001_5ESN and ALLe010_CR) the physical extent of the CGM is approximately 
a constant fraction of r s at all z. 


4.5 Exponential form of the column density profile and 
scaling with ionization energy of ions 

The column density profiles plotted on the log-log scale in 
Figure [6] exhibit a sharp turnover at a particular impact parameter. 
Although often interpreted as indication of an “edge” to the ion 
distribution, such sharp turnover can also simply reflect the smooth 
exponential decrease of column density. Indeed, the same profiles 
of the ALL_Efb_e0 01_5ESN run shown on the log-linear scale in 
Figure[7]can clearly be well approximated by an exponential profile 
(i.e., a straight line in the log-linear plot). Some ions exhibit two 


exponential components, but the outer component is at the column 
densities that are well below the sensitivity of current observations. 

The exponential profile N(d) = N 0 exp(— d/h s ), or equiv¬ 
alently log 10 N(d) = log 10 N 0 - hs ^ loy is characterized by a 
single scale, which we will call the scale height , h s , to differentiate 
from the exponential scale lengths of galaxy light distribution. 

The scale height of the profile provides a simple estimate of 
the extent of the CGM. |Nielsen et al.| ( |2013| ) and |Borthakur et al.| 
( 2015) show that the observed column density profiles of some of 
the ions can indeed be described by an exponential profile. Note, 
however, that in both studies the exponential distributions are fit- 
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Figure 6. Column density profiles for four representative low to high ions for our models and their evolution from z = 3.0 to z = 1.0. At high redshift 
(z = 2 — 3), the CGM seems to be insensitive to the feedback recipes except for the very weak star formation efficiency eff = 1%. At z = 1, the fiducial 
model starts to deviate from 5ffsNii and cosmic ray models. The extent of the profiles in units of R v i r shrinks with redshifts mainly due to the fast-evolving 
R v ir relative to the physical extent. 
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Figure 7. Column density profiles from low (HI and Mgll) to high (CIV and O VI) ions in the ALL_Efb_e001_5ESN run. Various lines are medians of the 
profiles at different redshifts (blue to red as low z to high z). Dark and light blue shaded regions represent 68% and 95% intervals for the profiles at z = 0, 
respectively. The left column shows the column density profiles in the inner regions remained relatively unchanged over 11 billion years at z ^ 3. This is 
why scaling with slower evolving scale radius, r s (right panel) brings profiles at different z into agreement better than scaling with R vir (middle column). At 
the same time, scaling with r s also accounts for the different masses of the galaxies, which is necessary to put a population of observed galaxies in an equal 
footing. Note that higher ionization ions have progressively more extended profiles. 
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Figure 8. Evolution of the virial radius and halo scale radius, r s . Virial 
radius is measured from the dark matter halo at each snapshot of the 
simulation ALL_E fb_e 0 01 _5ESN. The corresponding scale radius r s = 
^vir/cvir is calculated via the M v - ir — c v i r model in Diemer & Kravtsov] 
( 20151 . At high redshifts, the inner halo characterized by 4 r s tracks the 
growth of the halo virial radius, f? v ir, until z ~ 2. Afterwards, R vir con¬ 
tinues to grow, while the inner region of the dark matter halo approaches a 
constant value. 


Figure 9. Scale height of the exponential profile of the column den¬ 
sity as a function of the impact parameter predicted in our run 
ALL_Efb_e001_5ESN at z = 0 as a function of ionization energy of the 
corresponding ion. The left axis shows the scale height normalized to the 
halo scale radius, while the right axis shows the scale height normalized to 
the halo virial radius. Color points show results of the fits of the column 
density profile of a particular ion, while the dashed lines show the best fit 
power law relations (see eqs. |8|9) . 


ted to the equivalent width profiles instead of the column density. 
Due to non-linearity of the conversion between N and W r , scale 
heights from our fits should not be directly compared with ones 
from W r profiles. Instead, a conversion from N to W r is needed 
before comparisons (see Figure |TT|). 

The best fit scale heights of the run ALL_E f b_e 0 01 _5E SN for 
the inner exponential profiles of different ions are shown in Fig. [9] 
and presented in Table[2] Interestingly, Figure[9]shows that the scale 
height exhibits a tight power law scaling withion ionization energy: 
h s oc E™ n . Previous theoretical (e.g., |Hummels et al.|2013||Ford 
et al. 2015| ) and observational ( [Liang & Chen|20 14|| Johnson et al. 
2015 studies have also found that higher ions have more extended 
distributions compared to the low ions. Here we present the first 
quantitative characterization of this trend and demonstrate that it is 
a remarkably tight function of the ionization energy of ions. 

The best fit power law relation between scale height normal¬ 
ization to the virial radius, R v i r , or scale radius, r s , is characterized 
by the slope m and offset b : log 10 (x) = m log 10 Ei on + b, where 

for x — h s /r s : m = 0.742 ± 0.005, b Vs = -1.29 ± 0.02. ( 8 ) 

for x — h s /Rv ir : m = 0.742 db 0.005, 5# vir = b Vs — logc v i r 

(9) 

where c v ir ~ 10.5 is the concentration parameter of the halo at 
z = 0, computed from the M v i r — c v ir relation ([Diemer & Kravtsov| 
poT5) . We also find a similar power law slope m = 0.72 for the 
ALL_el0_CR at z = 1. We estimated the best-fit scale heights h s 
for each ion by a simple % 2 minimization fit to the simulated scaled 
profiles log Ni OU — d/r s . The fit is done to the range of radii, 0.1 < 
d/r s < 4, and column densities log N c < logiV < 17.2 cm -2 (21 
for HI), where log N c is the column density where we can clearly 
see the transition between the inner and outer exponential distribu¬ 


tions. The limits are chosen to exclude the ISM of the galaxy and 
limit the fit only to the inner exponential component. 

The errors on the parameters are estimated by bootstrapping 
the column density measurements of many lines of sight and thus 
account for the scatter in the profiles. We also varied the selected 
ranges of the data for the fit and find that h s can shift systematically 
by rsj 5 — 10%. Nevertheless, the power law dependence of h s 
on Ei on remains the same. We therefore caution readers that when 
comparing scale heights with observations, the data should be fitted 
in the same range. 

The trend of the scale height with ionization energy is likely 
due to two separate reasons. First, abundance of low ions is deter¬ 
mined by photoionization equilibrium and their distribution traces 
a particular range of the ionization factor, U, and thus a particular 
range of gas density profiles (see, e.g., §16.9 in [Mo et aL|2010) . 
Ions with higher ionization energies, such as CIV and O VI, can be 
created both by photo- and collisional ionization. Overall, they fa¬ 
vor regions with higher U and thus lower gas densities, which can 
be partly responsible for the trend we observe. The dependence of 
scale-height on the ionization energy of ions may thus be related 
to the fact outflowing plumes of gas expand as they move out to 
regions of lower density and pressure. 

To test this idea, we carried out an analysis with and with¬ 
out an ionizing background to compare the relative contribution 
from photo- and collisional ionization. Indeed, with the addition 
of photo-ionization, low ions are lifted into higher states in the low 
density regime (at large distances). This effectively boosts the scale 
heights of high ions, while decreasing the scale height of the low 
ions. Thus, the ionization factor dependence on density is at least 
partially responsible for the h s — Ei on scaling. 

The fact that higher ions arise from shocked gas that was car¬ 
ried out by older winds (e.g., [Ford et al. 12015) may also contribute 
to the correlation. An NFW-like potential is quite shallow in terms 


h s /R 
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Table 2. Scale Heights From Simulations 


Ionization State 

£ion [eV] a 

hs/ Rv\r 

h s /r s 

HI 

13.60 

0.042 db 0.004 

0.44 ± 0.04 

Mg II 

15.04 

0.025 db 0.001 

0.26 ± 0.01 

NV 

97.89 

0.142 ±0.012 

1.50 ±0.13 

Ne VIII 

239.09 

0.447 ± 0.006 

4.69 ± 0.07 

CII 

24.38 

0.047 ± 0.003 

0.50 ± 0.03 

cm 

47.89 

0.068 ± 0.003 

0.72 ± 0.04 

CIV 

64.49 

0.073 ± 0.004 

0.76 ± 0.04 

Sill 

16.35 

0.033 ± 0.002 

0.35 ± 0.02 

Si III 

33.49 

0.046 ± 0.002 

0.48 ± 0.02 

Si IV 

45.14 

0.053 ± 0.003 

0.56 ± 0.03 

OI 

13.62 

0.032 ± 0.002 

0.34 ± 0.02 

OVI 

138.12 

0.315 ±0.005 

3.31 ±0.05 

OVII 

739.31 

0.336 ± 0.004 

3.53 ± 0.04 

OVIII 

871.38 

0.659 ± 0.012 

6.91 ±0.13 


Morton 


a Ionization energy from a compilation of atomic data in 
Scale heights h s from exponential fits, log 10 7Vi on = log 10 JSIq — d/(h s ln(10)) 


(2003 


of the first component of the column density profiles for ALL_Efb_e001_5ESN 
at z = 0. We see a power law dependence h s oc fdP^ 4 . Note that 
h s /r s = c v ir x h s /R v i r , where c v i r ~ 10-5 is the concentration 
parameter, computed via M v - ir — c v i r model in 


Diemer & Kravtsov 


Uncertainties on h s are computed via bootstrap sampling. 


(20151. 



of its logarithmic slope out to r « 1 r s and steepens at large radii. 
It is thus tempting to associate the characteristic scale height of 
the low ions (« 2 r s ) with the region of the shallow potential out to 
which recent winds propagate easily and then stop when the deriva¬ 
tive of the potential steepens. However, we have checked the evolu¬ 
tion of gas outflows in our simulations and found that even the re¬ 
cent outflows in the ALL_Efb_e001_5ESN run propagate to con¬ 
siderably larger radii than r s . The scale height of low ions is thus 
likely determined by the overall shape of the gaseous halo density 
profile and the corresponding profile of the ionization factor. 

The prediction of the exponential profile and the tight correla¬ 
tion of the scale height with E- lori can be tested against observations. 
The correct treatment of the profile fitting in the observational data 
requires proper consideration of the upper limits and is out of the 
scope of our paper, but we plan to return to it in a future study. 


8 9 10 11 12 13 14 

log M vir 

Figure 10. Distribution of halo mass M v i r probed by various ions in the 
compiled low redshift data set. See references in Table[3] 


ing our adopted stellar mass-halo mass relation (see below). Fig¬ 
ure [ 10 ] shows that although the overall range of halo masses is 
broad, individual ions may be probed by a narrower range of halo 
masses. 

In this section, we first discuss details of how the comparison 
between observations and simulations is made. We then compare 
model predictions with observed column density, equivalent width 
and covering fraction profiles around galaxies and discuss the im¬ 
plications for galaxy formation and feedback. 


5 RADIAL COLUMN DENSITY DISTRIBUTIONS: 
SIMULATIONS VS. OBSERVATIONS 


5.1 Rescaling observations and simulations for fair 
comparison 


To compare our simulations with observations, we have com¬ 
piled a set of existing measurements of the CGM absorbers probed 
with different ions around galaxies with stellar mass estimates. The 
data set is summarized in Table [3] which specifies the number of 
galaxies, their mean redshifts and stellar masses, and the range 
of impact parameters probed by a particular sample. The absorber 
galaxy samples span a wide range of stellar masses from dwarf 
galaxies with M* « 10 7 M© to L* galaxies with M* « 10 11 M© 
at z ~ 0, and Lyman-break galaxies with (M*) ~ 1O 11 M 0 at 
z ~ 2. Furthermore, Figure [lO] shows the ranges of virial halo 
masses for the galaxies which generated a particular ion absorber, 
where virial masses were estimated from galaxy stellar masses us- 


Observations of low-z absorbers in our data set span more than 
four decades of stellar mass ( Chen et al . 2010, Tumlinson et al.| 


2011||Liang & Chen|2014||Werk et al.|2014| Bordoloi et al.|2014[ 

Johnson et al. 2015), while our simulations model a single progen¬ 
itor of mass close to L* galaxies. We thus need to rescale both data 
and simulations in a way that makes the comparison sensible. 

In the previous section, we showed that the column density 
profiles in simulations scale most optimally across redshifts not 
with the virial radius of their parent halo, but with the scale radius 
of their halo density profile, defined as the radius where the pro¬ 
file has logarithmic slope of —2. To test whether a similar scaling 
applies to observed profiles in Figure |TT] we show the equivalent 
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Table 3. Observation data sets of CGM 



# of Galaxies 

<*> a 

Impact Parameter Range (kpc) 

log M*/M e b 

Transitions 

Johnson et al. 2015 

71 

0.271 db 0.088 

63 - 991 

10.14 ±0.92 

Lyo:, O VI 

Liang & Chen 2014 

195 

0.041 db 0.044 

32 - 499 

9.92 ± 1.20 

Lya, CII, CIV, Si II, Si III, Si IV 

Bordoloi et al. 2014 

43 

0.027 zb 0.023 

14 - 135 

9.5 ± 0.48 

CIV 

Werk et al. 2013 

44 

0.221 db 0.053 

18 - 154 

10.61 ±0.50 

Lya, CII, Sill, Si III, Si IV 

Tumlinson et al. 2011 





OVI 

Chen et al. 2010 

75 

0.239 db 0.094 

8.5 - 119 

9.9 ± 0.58 

Mg II 

Steidel et al. 2010 

512 (stacked) 0 

2.2 db 0.3 

10 - 125(280) d 

9.85 ± 0.46 e 

Lya, CII, CIV, Si II, Si IV 


“Median redshift and dispersion of the galaxy samples. 

^Median stellar mass and dispersion of the galaxy samples. 

“Mean equivalent width W r measured from stacked background galaxy spectra. 
^Maximum impact parameter is 280 kpc f or Lyo: and 12 5 kpc otherwise. 


“Median and dispersion of stellar mass for Steidel et al. (2010», a representative sample of galaxies in Reddy et al. (2012). 



d [kpc] df Ry ir dfv s 

Figure 11. Equivalent width distribution with impact parameter for the low- 2 : galaxies observed by the COS instrument onboard HST (see Table [3] for 
description of observational samples and references) and z ~ 2 stacked measurements by [Steidel et al.|(2010| i. Blue pentagons are detections, while gray 
pentagons with arrows are upper limits. The color 2D histogram shows distribution of equivalent widths for the lines of sight through the halo of the MW-sized 
galaxy in the ALL_Ef b_e 0 01 _5ESN run, assuming turbulent broadening of 6nt = 20 km s _1 . Black dashed lines represent W r profiles converted from the 
best fit exponential profiles to the simulated column density profiles for corresponding ions. The left panel shows profiles against the impact parameter d in 
physical kpc, while in the middle and right panels d is rescaled by the virial radius, R v i r , and halo scale radius, r s , respectively. 
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Figure 12. Equivalent width profiles for three representative ions in the run with cosmic ray feedback at z = 1 as a function of impact parameter scaled by 
the scale radius, d/r s (2D color histogram) compared to observations (points; the point types are the same as in FigurefTT). 


width profiles of absorbers in the ALL_Efb_e0 01_5ESN simula¬ 
tions along with observational data at redshifts z £ 0 — 0.3 and 
z ~ 2 with radii in physical kpc, and rescaled by the virial ra¬ 
dius of host halo R v i r (z) and r s (z). We show evolution of the W r 
profiles for CIV and Si II only, because the profiles for these ions 
exhibit some of the sharpest “edges” and because results for other 
ions are similar. 

We convert our predicted column density into rest frame 
equivalent width W r using the curve of growth analysis discussed 
in section 3.3. To account for turbulence not fully resolved in simu¬ 
lations in the relatively low resolution regions of gaseous halos, we 
use the b parameter as b * 2 = b% + b ^ T , and adopt 6 nt ~ 20 km s -1 
(see also |Oppenheimer et al.|2012| , although we checked that using 
&nt also provides a reasonable match to observations. 

In the simulations the virial radius of the galaxy is known from 
the density profile of the matter distribution. In observations, how¬ 
ever, we need to estimate it statistically using stellar masses. We 
do this using the average stellar mass-halo mass relation derived by 
|Kravtsov et al. 1(2014) using the abundance matching. This relation 
uses the measurement of galaxy stellar mass function by Bemardi 
|eTaTlp0T3] ), which is based on updated photometry of the SDSS 
main galaxy sample that corrects for background subtraction errors 
in the previous SDSS data releases. 

We assume a non-evolving M* — M v i r (z) relation (e.g., 
|Behroozi et al.|[20l3| to convert M* to M v i r for each observed 
galaxy in the dataset. We then use M v i r to compute virial radius us¬ 
ing its definition itVr(z) = [3M v ir(^)/(47rA v irp(z))] 1//3 , where 
Avir(z) is the redshift-dependent “virial” density contrast com¬ 
puted using approximation to the spherical collapse overdensity 
given by Bryan & Norman (1998 ) jThe scale radius is obtained 


from the virial radius assuming the median concentration of halos 
of mass M v ir at redshift z, predicted by the model of Diemer & 
Kravtsov <|2015|): r s — #vir/c v ir (M v i r , z). 


Figure [TT] shows that evolution in the physical radius from 
z & 2 to the present epoch in the observational sample is quite mild 
(as found previously by |Chen|2012||Liang & Chen|2014| ). The fig- 


2 Adopting other commonly used definitions of the “virial radius”, such as 
#200m or #2000 would not qualitatively affect our conclusions. 


ure also shows that the drop in the W r profiles is sharper when r s 
rescaling is used instead of # v ir- Furthermore, the profile rescaled 
by r s does not evolve significantly between z ~ 0 and z ~ 2. 

This alignment of the profiles at different z is remarkable be¬ 
cause properties of the low -2 and z ~ 2 populations of galaxies are 
drastically different. For example, typical star formation rates dif¬ 
fer by over an order of magnitude. Given the theoretical arguments 
for r s scaling and the fact that it appears to work for real CGM, we 
will scale impact parameters in the simulations and observational 
sample by the radius of 4 r s in the comparisons presented below. 
We choose the multiple of four because at z ~ 2 this radius is actu¬ 
ally close to the virial radius of galaxy halos (see |More et al.|2015| 
for a detailed discussion). 


5.2 Comparison of the column density profiles in simulations 
and observations 

Figures [13] [14] and [15] compare column density distributions 
with d/Ar s in the two simulations that were run to z — 0, 
ALL_Efb_eO 10 and ALL_Efb_e001_5ESN, and observations for 
all commonly observed transitions. It is immediately apparent that 
the absorber column densities in the fiducial run ALL_Ef b_eO 10 
with eff # 10% for all ions are consistently low by up to sev¬ 
eral orders of magnitude compared to observations. Interestingly, 
Agertz & Kravtsov (2015 2016 ) show that this run predicts prop¬ 
erties of the stellar and cold HI gas component of the galaxy in 
very good agreement with properties of a typical late type galaxy 
of this mass. Thus, there is a striking difference in how successful 
this run is in producing realistic stellar components of galaxies and 
the discrepant gaseous CGM it predicts. 

In contrast, the ALL_E f b_e 0 01 _5E SN run matches the CGM 
column density distribution of all ion species in the log N — d/Ar s 
plane reasonably well, especially given that we are comparing a 
single object to galaxies of a range of stellar masses. However, 
this run produced stellar component of the galaxy with a spheroidal 
morphology, exponential stellar surface density profile and unreal¬ 
istically large half-mass radius. This shows that the CGM provides 
a stringent and orthogonal constraint on the galaxy formation sim¬ 
ulations relative to the optical observations of galaxies. 

It is interesting that the CGM in the ALL_Efb_eO01_5ESN 
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Figure 13. Comparison between predicted radial column density distributions of various ion absorbers (2D histograms) with observations at z ~ 0 — 0.3 
(pentagon points; see Table [3j. Blue points show detections, while gray points show non-detections with downward arrows denoting 2cr upper limits. The 
ion for which the profile is constructed is indicated in the upper right. The left column shows profiles predicted by the fiducial model with all feedback + 
-E/fb, eff = 10%. The right column shows the predicted profiles from the model with all feedback + = 1%, 5-Esnii- The 2D histograms show that 

the majority of the points in ALL_Efb_eO 10 are many orders of magnitude below observations, while profiles in the ALL_Efb_e001_5ESN run are closer to 
observations. 
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Figure 14. Similar to Figure[l3] but for Si II, Si III and Si IV. 
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Figure 15. Similar to Figure[l3] but for Mg II, O VI and Ne VIII. 
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run reproduces the sharp drop in the column density of absorbers 
apparent for all transitions, which coincides with the location where 
a large fraction of measurements are upper limits. [Liang & Chen| 
(2014) have estimated this drop to be at d/R v i r ~ 0.7 for Sill, 
Si III Si IV, CII and CIV. A similar drop of the CIV column densi¬ 
ties was reported by Bordoloi et al. (2014) at d/R 2 oom ~ 0.5, and 
for O VI by Johnson et al. ( |2015| ) at d/R v i r ~ 1. 

The difference between d/i? 2 oom = 0.5 and d/R v i r = 0.7 
found for low and intermediate ions is simply due to the difference 
between R V1T (A « 360) and it^oom definition and the specific 
adopted M* — relation. Once consistent definitions are adopted, 
the two samples of data are in agreement. Indeed, if we re-compute 
R v i r using the Bryan & Norman approximation and |Kravtsov et al.| 
< |2014| ) M* — Mh for all galaxies compiled in this study, the low- 
ionization metal boundary is consistent with d/R v i r ~ 0.7. 

[Liang & Chen| j2014| found a hint that the CIV extends be¬ 
yond dj -R v ir = 0.7, but based only on absorbers around two galax¬ 
ies. Here, combining the data sets ofjLiang & Chen|(|2014]) and[Bor-| 
[doloi et al.| ( 201 ?), we see a stronger indication that CIV indeed 
extends to d « 6 — 7 r s ~ 0.9i? v ir and its distribution is thus more 
extended compared to the low-ions (e.g., Sill and MgII). Taking 
into account the result of | Johnson et al.| ( [2015| that the O VI distri¬ 
bution extends to d « i? v ir, observations indicate that the extent of 
the absorbing gas increases with increasing ionization energy of the 
ions. This is qualitatively consistent with the tight trend between the 
scale height of the exponential column density profile and ioniza¬ 
tion energy of the corresponding ion that we found in simulations 
(see Figure [9] and Table [2]). It remains to be seen whether observa¬ 
tions will confirm the tight correlation of the profile scale height 
with ionization energy predicted by our simulations. 


form a bootstrap re-sampling of the simulations to estimate the un¬ 
certainty of the covering fraction profile in simulated galaxies. 

Figure |T6] shows the covering fraction profile predictions for 
MgII, CIV and OVI for runs ALL_Efb_e010, ALL_e010_CR 
and ALL_Efb_e001_5ESN. As discussed above, the column den¬ 
sities of most lines of sights in the ALL_Efb_e010 run is orders 
of magnitude too low compared to observations. These low column 
densities are below the detection limit of the observations, resulting 
in low covering fraction. This is apparent for Mg II, CIV and even 
OVI. The only region with non-zero k is with distance d < r s 
(or ^ O.litLir for the simulated halos), at the extended disk where 
density is high and temperature is low. 

On the other hand, the ALL_Efb_e001_5ESN run produces 
cold/warm outflows that reach to larger distances, producing a 
higher covering fraction. The covering fraction profile predicted in 
this run is much closer to the profile derived from observations. 
Although the match is not perfect, we note that we are compar¬ 
ing a single object in simulations to a sample of galaxies spanning 
a wide range of stellar masses. Thus, meaningful conclusions can 
be drawn only from larger samples of simulated galaxies. Never¬ 
theless, the fact that simulations reproduce high covering fraction 
extent to large radii in rough agreement with observations is en¬ 
couraging. 

The run ALL_e010_CR, which includes cosmic ray feedback, 
predicts profiles at z = 1 (the lowest z to which the simulation 
was run) that are intermediate between those of ALL_Ef b_e010 
and ALL_Efb_e001_5ESN. For OVI it is considerably closer 
to ALL_Efb_e001_5ESN. This indicates that cosmic ray-aided 
winds can help to bring predictions in better agreement with ob¬ 
servations. 


5.3 Comparing the CGM Covering Fractions 

The scatter of the column density in the region of sharp drop of 
ion column densities shown in Figures[l3][T4]and[T5]is large, and at 
some radii may span orders of magnitude. The large scatter means 
that it may be misleading to compare median or average profiles 
and one has to compare distributions of column densities directly, 
as we did above. Another way to characterize the column density 
distribution is via the covering fraction k, which measures a frac¬ 
tion of area covered by absorbers with column density or equivalent 
width larger than a given threshold. 

Different detection thresholds reported for different observa¬ 
tions present a challenge for a uniform comparison. Some stud¬ 
ies adopt rest frame equivalent width of Wo = 0. 05 — 0.1 A 
and log No — 13.5 cm -2 as detection thresholds |Liang & Chen 
|2014 , Johnson et al. 2015). We adopt these values here for compar¬ 
ison with the derived covering fraction of data from Liang & Chen 
< |2014| ) and | Johnson et al.| ( |2015| ) for OVI. For observations where 
only e quivalent width measurements are available, Liang & Chen 
42014 4 adopted the threshold W r = 0.05 A. We convert this value 
from Wo to log Ao = 13.75 for the calculation of the covering 
fraction, using their estimated b parameter of 29 km s -1 . 

In principle, one can use a different way of treating the up¬ 
per limits in which the upper limit is considered either as detection 
with the corresponding column density or equivalent width, or as 
effectively zero column density. The former choice gives the largest 
possible covering fraction, while the latter gives the smallest pos¬ 
sible covering fraction. These two limits are shown as blue bands 
in the panels of Figure |T6] Overall, the two methods give estimates 
consistent with each other, but the second approach in some cases 
results in a wider range of allowed covering fractions. We also per- 


6 DISCUSSION 

In the preceding section we have presented comparisons of 
properties of gaseous halos forming around a progenitor of a « L* 
galaxy in simulations with different parameters and implementa¬ 
tions of the star formation - feedback loop. We find that properties 
of the gaseous CGM, in particular, spatial distribution of various 
ions, are highly sensitive to the details of these implementations and 
parameter choices. For example, the fiducial simulation of |Agertz| 
& Kravtsov (2015 2016 ), which produces a very realistic central 
galaxy with the correct stellar mass, size, angular momentum, rota¬ 
tion curve, bulge-to-disk ratio, stellar and gas surface density pro¬ 
files does not produce extended CGM and is in striking discrepancy 
with observations. This illustrates that properties of galaxies and 
properties of their CGM provide strong complementary constraints 
on the processes governing galaxy formation. 

Variations of the stellar feedback model, such as adding feed¬ 
back modelling due to cosmic rays or simply changing parameters 
governing star formation and stellar feedback, affect properties of 
gaseous halos around simulated galaxies and produce a more ex¬ 
tended CGM in better agreement with observations (see, e.g., Figs. 
|13|15| above). This is particularly true for ions with higher ioniza¬ 
tion energies, such as C VI and O VI, indicating that these simula¬ 
tions correctly capture the thermal and density structure and metal- 
licity profile of hotter gas. 

There are indications that the distribution of low-ions in all of 
our simulations is somewhat less extended than indicated by obser¬ 
vations. In the fiducial run, this is due to the lack of low tempera¬ 
ture gas at large radii at low z (see also |Hummels et al.|2013| ). This 
is mainly because in the fiducial run there are no global winds at 
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Figu re 16. Comparis on of the predicte d covering fraction profiles with measuremen ts of Mg II A2796 (Chen et al.| 20 To} , CIV A1548 ( [Liang & Chen|20l4l 
and O VI A1031 ( Tumlinson et al.| 2 oTT Johnson et al.|2015 9 transitions. Detection threshold Wo (CIV A1548) = 0.05 A is adopted 
< |2014 1 . We convert the threshold to column density using a cur ve-of-growth analysis adopting the doppler parameter b = y/2a = 


B ordoloi et jdj2014 


Liang & Chen 


Johnson et al. 


( |2015jf . log A 0 (MgIIA2796) = 13 cm - 2 is chosen 


from 

29 km s _1 . Detection threshold log No(0 VIA1031) = 13.5cm -2 is adopted from|J _ 

based on the level upper limits in |Chen et al. 1(2010) . Error band represents 68 % confidence level bracketed by two limiting cases of consideration of upper 
limits above the detection threshold (e.g., as detection and as non-detection). In all three ions, ALL_Efb_e001_5ESN and ALL_e010_CR predictions agree 
with the inner region (< 2 — 4 r s ) of the data while under-predicting the covering fraction of the outskirts of the halo (~ 10r s ). Note that ALL_e010_CR is at 
z = 1 while the other runs are at z = 0. It is conceivable that halo in the CR run will continue to grow at lower redshift, as suggested by Fig. [5] 


z < 0.5 — 1 and thus there is no cold/warm gas lifted from the disk 
of the galaxy into the halo. In contrast, in the 5Esn run, the winds 
continue to z = 0 and therefore the CGM is filled with cold/warm 
gas at low redshifts. 


The low metallicity of the halo in that run may also play a role 
in the low column density of absorbers. For example, the total col¬ 
umn density of all carbon is still lower than the observed column 
density of CIV alone by a factor of ~10. We think this is because 
the winds at high z in the fiducial run are so strong that they dis¬ 
perse metals over a large volume, lowering the overall metallicity. 


Some recent observations indicate the presence of low-ion ab¬ 
sorbing clouds with sizes below the resolution of the simulations 
(Muzahid 2014, Crighton et al. 2015), which is particularly poor in 
the gaseous halo far away from the dense regions of the disk. It is 
possible that a sub-grid model of small clouds created by nonlinear 
thermal instabilities ( [Joung et al.[2012|) can ease the disagreement. 
Gas inflow along filaments and feedback-driven outflows could 
plausibly excite such instabilities, which then could create clouds 
that would boost the amount of cold/warm gas present. Some low- 
redshift observations, on the other hand, indicate large physical 
sizes (Davi s et ak]|2015|) and non-hydrostatic states of many ab¬ 
sorbers ( Werk et al.|2014| consistent with the properties of recent 
outflows of the kind we observe in our ALL_Efb_e0 01_5ESN and 
ALL_e010_CR runs. Low-ion absorbers in galaxies thus may orig¬ 
inate both from large-scale outflows and from small clouds forming 
by thermal instability that outflows help to excite. 


Overall, the comparison of simulation results with different 
star formation and feedback models has produced valuable insights 
about the radial distribution of absorbers, which we discuss below. 
Below we will also discuss how the results of our simulations com¬ 
pare to other recent theoretical studies on this subject. 


6.1 Column density profiles of the CGM absorbers: is there 
an ion boundary of the CGM? 

As we discussed in section H3 observations of absorption 
lines of a variety of ions, such as Si II, Si III Si IV, CII, CIV, and 
O VI, exhibit a sharp drop in the incidence of absorber detections 
and their covering fraction beyond a certain radius ( [Liang & Chen| 
|2014[|Bordoloi et aL|2014[|Johnson et al.|2015| ). Some researchers 
have interpreted this drop-off as an “ion boundary.” 

Our simulations of a ~ L* galaxy progenitor at different red- 
shifts predict approximately exponential column density profiles 
for all ions. When plotted on a log-log scale, the exponential pro¬ 
files exhibit a sharp turnover at approximately the “half-mass ra¬ 
dius” or « 1.68 h s , where h s is the scale height of the exponen¬ 
tial profile. Thus, our simulations indicate that the perceived sharp 
turnover in the incidence of high column density absorbers is sim¬ 
ply a manifestation of a continuous underlying exponential profile. 

It is worth noting that de-projection of the projected exponen¬ 
tial column density profile via the Abel integral gives the modified 
Bessel function of the second kind, K 0 (x) for the radial profile of 
3D number density (see, e.g., Appendix in |Pitts & Tayler||1997| ). 
The Ko(r) profile has a shallow slope at small r, which steepens 
at larger r and we find that the K 0 (x) function does approximate 
the gas density profile in our simulations at r ~ 20 — 200 kpc 
quite well. The origin of the approximately exponential total col¬ 
umn density profile, Nn(d), is thus in the shape of the 3D gas den¬ 
sity profile. 

The column density profiles of specific ions are additionally 
shaped by the metallicity profile and by the ionization factor and 
thermal structure of the gas that affect their photo- and collisional 
ionization rates. These factors are ion-specific and give rise to pro¬ 
files that vary systematically with the ionization energy of ions. 
The metallicity profile in the central regions varies in different sim¬ 
ulations. In the ALL_Efb_e0 01_5ESN run the metallicity rapidly 
increases with decreasing radius at r < 50 kpc, while in the fidu¬ 
cial ALL_E f b_e 010 run the metallicity remains low and its profile 
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Figure 17. Column density profiles of Mg II, CIV and O VI in the outer 
region of the hot gaseous halo share the same underlying exponential profile 
traced by the total gas column density profile. All other ions (not shown) 
also show similar exponential profiles. Column densities for ions have been 
offset by log Nq = 6. 


remains flat down to small radii. The difference is that the former 
simulation drives outflows at low redshifts that enrich the central 
region 50 kpc) and provides a steady supply of fresh warm gas 
that gives rise to low-ion absorbers. 

We find that column density profiles of all ions trace the over¬ 
all shape of Nu (d) quite closely at d > 50 kpc, as shown in Fig¬ 
ure [IT] because the metallicity profile at these radii in our simula¬ 
tions is quite flat. The low-ions at the large radii have column den¬ 
sities much below the sensitivity limits of observations. However, 
higher ionization ions, such as O VI, do produce detectable absorp¬ 
tion. Thus, the exponential column density profile of ions like O VI 
simply reflects the underlying total gas profile. At r < 2 r s ~ 50 
kpc the profiles are more difficult to interpret, but it appears that the 
cycling outflows at these radii do result in a metallicity profile that 
is also approximately exponential in projected radius. Other fac¬ 
tors, such as dependence of the U -factor on radius, probably also 
contribute to the shape of the column density profiles of specific 
ions. In particular, the U-factor almost centrainly plays a role in 
setting the tight correlation of scale height with ionization energy 
we find for the low ions. 

Encouragingly, there have been a few recent studies that show 
exponential profiles are a good fit to the observed data. Study by 
|Nielsen et al.[ ( [2013) ) finds that an exponential profile provides a 
better fit to the Mg II equivalent width profile, as compared to a 
power law. Borthakur et al. (2015) finds an approximately expo¬ 
nential form for the Lyo; equivalent width profile. Bor doloi et al.| 
< |2014| ) have also used an exponential function to describe the sharp 
drop in the equivalent width of CIV absorber column density in 
their sample, although overall they used a more complicated pro¬ 
file given by a product of a power law and an exponential function. 


6.2 Column density profiles of the CGM absorbers: scaling 
with halo mass and redshift 

Interestingly, as shown in section[5]above, the collection of ob¬ 
servations considered in our study indicates that the scale height of 
the column density profiles for observed absorbers around galaxies 
of a wide range of stellar masses occurs at approximately the same 
multiple of the halo scale radius, r s at different redshifts. The tran¬ 
sition in galaxies at z ~ 2 in the sample of |Steidel et al.| ( [ 20 T 0 | also 
occurs at the same multiple of r s . Regardless of the precise form 
of the profile, this scaling is consistent with the scaling predicted in 
our simulations for the inner regions of the CGM profiles. 

The scaling of the overall column density profile with r s is 
due to the slow evolution of the inner gaseous profiles. Simula¬ 
tions of halo evolution in CDM models show that the inner pro¬ 
files of galaxy-sized halos are in equilibrium and evolve little after 
1 — 2 (Prada et al.|2006[ Diemand et al.|2007]|Cuesta et al.| 


2008||Diemer et alj 2013, Zemp 2014, Diemer & Krav tsov|2014| 


Correa et al. 2015), while the profiles in the outskirts (r > ifeooc) 


continue to accrete mass and evolve ( [Diemer & Kravtsov||2014[ 
More et al. 2013] see also discussion in Section |4~4| above). The 
virial radius of the halo tracks the evolution of the outer profile and 
the outer splashback boundary of the halo fairly closely and thus 
evolves considerably faster than the profile in the inner regions. 
The evolution of the latter can be better characterized by the halo 
scale radius ( [More et al.|2015| ), which also exhibits slow evolution 
at low redshifts ( [Bullock et al.|2001~| ). If the overall gaseous profile 
is in approximate equilibrium, this can explain the apparent scaling 
of the column density profile evolution with r s . 

Note that at a fixed redshift , the scalings with R vir and r s 
are nearly equivalent because r s = Rvir/c v i r and concentration 
is a very weak function of halo mass: c v ir oc oc Rf®' 3s , 

where s ^ 0.08 — 0.1 at z = 0 and is shallower at higher redshifts 
(e.g., |Bullock et al.||200l| |Neto et aL]|2007[ |Diemer & Kravtsov] 

|2015| >. Thus scaling with r s at a fixed z is equivalent to scaling 
with R^- 1 ' 25 . At the same time, concentration exhibits a rather 
strong evolution with redshift due to fast evolution of R v i r and thus 
redshift evolution of r s and R v i r and corresponding scaling as a 
function of redshift will be very different. Therefore, our results 
indicate that for samples spanning a wide range of redshift, scale 
radius r s should be used to rescale radial scales, such as the impact 
parameter. 

The remaining question is why the column density profiles 
shaped by stellar feedback driven outflows in simulated galaxies 
scale with r s . We believe this can be understood from the following 
considerations. The specific energy of the wind when it is ejected is 
E = v^/2 + 0(r w ), where v w is initial wind velocity at the launch 
radius, r w and is gravitational potential at r w . If the outflowing 
gas is gravitationally bound (E < 0), the wind will stop at some 
radius r ou t < Rv ir- If we assume that energy losses due to ram 
pressure force from the tenuous halo gas on the wind can be ne¬ 
glected, then the wind stops at E « fi(r ou t). If r ou t is sufficiently 
large we can approximate the potential by the form expected from 
the NFW profile: 


Hr out) » — 4.63 V r 


2 

max 


ln(l + Xout) 


%out 


( 10 ) 


where V m ax is the maximum circular velocity of the halo, x ou t = 
r out / r s and r s is the scale radius of the NFW profile. Although the 
potential in the region where stars launch winds is likely not de¬ 
scribed by the NFW form due to contributions of the baryons in the 
galaxy, its overall amplitude should still scale as oc — V^ ax . 

At the same time, for both momentum- and energy-driven winds 
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we expect v w oc V ma x ( [Murray et al.|2005||Booth et al.|2013[|Mu-| 
|ratov et al.|2015|>, the scali ng that is also indicated by observations 
( [Schwartz & Martin|2004||R~upke et al.|2005| [Schwartz et al.|2006[ 

although see Heckman et al. 2015, submitted). Therefore, the above 
equations suggest that 

—-- « Cl + c 2 ^(r 0 ut), (11) 

%out 

where a , C 2 are constants and ^ (r ou t ) does not depend explicitly 
on Vmax. This equation means that x out is independent of V ma x and 
hence halo mass and thus r ou t oc r s . Although individual fountains 
may have a range of initial velocities and launch radii, from the 
above considerations we can see that the characteristic extent of 
the profile they shape may naturally scale with r s . 


6.3 The Effects of Self-Shielding, Non-Equilibrium and 
Local Starburst Radiation 

In this work, we use three simplifying assumptions in the cal¬ 
culation of ionization abundance: we assume that gas is optically 
thin, that it is in ionization equilibrium and ionized by the cosmic 
UVB only (i.e., not accounting for the local radiation from the stars 
in the simulated galaxy). We qualitatively discuss their effects and 
implications on the CGM here. 

First, self-shielding may increase the column density of low 
ions compared to our estimate. We find, however, that self-shielding 
can only change the column density significantly near the star¬ 
forming disk where density is high. For example, [Rahmati et al.| 
( |2013| ) shows that self-shielding affects gas at column densities 
greater than logiV « 17 cm -2 and our Figure [ 7 ] shows that this 
corresponds to r < 25 kpc. 

Second, the ionizing radiation from local star-forming regions 
will likely increase the ratio of low to high ions due to a softer 
stellar spectrum compared to QSOs. Assuming escape fraction of 
/esc = 1 — 2%, the contribution of local radiation is negligible at 
radii beyond « 30-50 kpc ( |Shen et al.|2012| >. Note that adding local 
radiation would affect our result in the opposite direction from self¬ 
shielding because it would increase ionizing radiation in the regions 
of high gas column density. 

Finally, deviations from photo-ionization equilibrium may 
boost ionization at lower temperatures. However, this effect be¬ 
comes less important at lower metallicity, especially in the pres¬ 
ence of photo-ionization by a radiation background ( |Oppenheimer| 
|& Sc haye 2013), the regime of our analysis. 

Given the current small observational statistics, and lack of 
strong arguments that these effects are significant at large radii, ne¬ 
glecting them should not affect our conclusions significantly. 


6.4 Comparisons with Previous Studies 


A number of recent studies have explored predictions for the 
observable properties of the CGM in cosmological simulations of 
galaxy formation. Some studies have focused exclusively on pre¬ 
dictions at high redshifts (z 


Voort et al. 

2012 

Shen et al.||2012 2013, Pallottini et ah 2014, 

Suresh et al. 

2015 

jFaucher-Giguere et al.|2015). Much of the high- 


z work has focused on exploring predictions for the distribution and 
observable properties of neutral hydrogen and dependence of pre¬ 
dictions on the implementation of feedback and associated winds 
(e.g.,|Barnes et al.|2011||van de Voort et al.|2012[|van de Voort &| 

|Schaye|2012[ Pallottini et al.|2014[|Faucher-Gigu£re et al.|2015| >. 
Studies that systematically analyzed properties of the CGM at 


low redshifts used two kinds of simulations. |Ford et al.| ( |2013[|2014[ 
2015) used a statistical sample of galaxies formed in a simulation 
of 32h~ 1 Mpc and 48 h~ x Mpc boxes, but which modelled the 
wind launching process phenomenologically by imposing a partic¬ 
ular specified wind scaling with halo properties (e.g., wind velocity 
and loading factor) as a function of halo mass. In these studies the 
authors have placed particular emphasis on the origin and dynam¬ 
ical cycles of the CGM. |Ford et al.| ( |2013[|2014] ) showed that low 
ions, such as Mg II, trace dense gas close to galaxies that were part 
of recent outflows and which will re-accrete onto a disk on a ~Gyr 
time scale, while high ions, such as O VI, have more extended dis¬ 
tributions and originate from ancient outflows (see also|Shen et al.| 
i2012). Our results are in qualitative agreement with the results of 
these studies. Moreover, the column density profiles presented in 
|Ford et al~] ( |2014| > and |Ford et al. 1(2015) have a form in qualitative 
agreement with results of our simulations. In particular, their pro¬ 
files are approximately exponential in the inner regions, with some 
ions exhibiting two exponential component profiles (see, e.g., Fig. 
10 in |Ford et al.|2013[ and Fig. 7 in |Ford et al.|2014| >. 

In the second type of simulations, which includes simula¬ 
tions presented in this work, formation of individual galaxies 
was modelled using different numerical codes (Enzo, RAMSES, 
Gasoline) in the “zoom-in” simulation where all mass and spa¬ 
tial resolution is focused on a single object ( Stinson et al.||2012| 


Shen et al.|2013||Hummels clal^ Ol 3 Ros karet al.|20141 |M arasco 


et al.|2015) . The higher resolution in the disk region in such sim¬ 
ulations allows a more sophisticated treatment of stellar feedback, 
which allows treating the wind launching process self-consistently 
without imposing particular wind scaling properties. Not surpris¬ 
ingly, the CGM properties in such simulations were shown to be 
quite sensitive to the details of star formation and stellar feed¬ 
back implementat ion ([Hummels et al.| [2013, Ro skar et al.| l2014, 
|Suresh etal.|20 15 1 Marasco et al. 20 15]. Again, this is in agreement 

with our results, which indicate strong dependence of the proper¬ 
ties of the CGM on the feedback parameters and implementation. 
We should note that the high sensitivity of the CGM properties to 
wind modelling is specific to the zoom-in simulations with self- 
consistent wind launching. [Ford et aT] ( |2014[ |2015| ) did not find 
significant differences in the CGM produced by the different phe¬ 
nomenological wind models used in their simulations. 

In Figure [18] we compare the CGM profiles from two of 
our simulations with SN feedback with the profiles of Hum- 
mels et al. (|2013) ^The figure shows that the profiles from our 


ALL_Efb_e010 are in good agreement with the “Medium Feed¬ 
back” model of |Hummels et ah ( 2013|), which produce a CGM 
inconsistent with observations. |Hummels et al.| ( |2013| ) found that 
artificially delaying cooling after gas is heated by SNe makes 
wind launching more efficient and produces a more extended CGM 
halo, although the extent of the ion distribution is still somewhat 
smaller than indicated by observations. In all transitions, our model 
ALL_Efb_e0 01_5ESN can be viewed as a similar variation of the 
feedback parameters, although instead of suppressing cooling, it 
assumes larger energy release per supernova. Physically, this can 
correspond to a larger fraction of thermal energy injected by a su¬ 
pernova retained by the gas or to top-heavy IMF. Our profiles for 
this simulation are indeed closest to the “high feedback” run of 
|Hummels et al.| ( [2013| ), although our profiles are somewhat more 
extended and are thus in better agreement with observations. 


3 The profiles are taken from the website of the authors: 

http://chummels.org/CGM.html 
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Figure 18. Comparisons of predicted column density profiles from low to high ions between our simulations with the profiles from simulations of |HummeIs| 
|et al.|p0l3| denoted as H+13). The profiles in our fiducial run are closest to the medium feedback run of H+13. 


One may question the physical validity of the cooling de¬ 
lay or increasing energy per supernova. However, we think that 
these choices may approximate some other physical processes 
that result in more efficient wind driving. For example, when we 
include cosmic rays with isotropic treatment of their diffusion, 
simulations produce CGM profiles quite close to those of the 
ALL_Efb_eO01_5ESN run, albeit without assuming larger energy 
per supernova. Overall, these results (and previous studies that pro¬ 
duce realistic CGM) show that feedback processes that contribute 
to wind launching should be quite efficient. 

As we discussed in the previous section, the column density 
profiles of the CGM ions are predicted to scale linearly with halo 
scale radius, r s , although we could naively expect that strong feed¬ 
back could break such self-similar scaling with halo parameters. 
Other recent studies have reported similar findings. For example, 
| van de Voort & Schaye| ( |2012|) show that outskirts of gaseous ha¬ 
los of galaxies in the OWLS simulation evolve as expected from 
the “virial” scaling relations. [Pallottini et al.| ( |2014] > find that CGM 
profiles around high-redshift galaxies (z = 4) are self-similar once 
scaled with itLir- They have found that the HI column density pro¬ 
files as a function of normalized impact parameter, d/R v i r , have 
power law form, while we find exponential form (or double ex¬ 


ponential) for all ions and neutral hydrogen absorbers. Our results 
indicate that column density profiles of the CGM ions do exhibit 
self-similarity, but they scale better with halo scale radius, r s , rather 
than with the virial radius. 


7 CONCLUSIONS 

We presented a detailed analysis of the observable properties 
of the circumgalactic medium in a suite of high-resolution cosmo¬ 
logical galaxy re-simulations of a Milky-Way sized halo with a va¬ 
riety of star-formation and feedback models. We have also con¬ 
trasted these predictions with a large set of existing observations of 
absorbers around galaxies of different mass from z ~ 0 to z « 2. 
Our specific findings can be summarized as follows. 

1. Our simulations indicate that the CGM probed by absorbers 
of different ions arise in gas that was ejected from the disk by stel¬ 
lar feedback. At low redshifts, the outflows in the MW-sized halo 
modelled in our simulations are not sufficiently energetic to leave 
the halo and instead produce plumes that reach a finite radius and 
then turn around to infall back onto the disk. We show that low 
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ionization absorbers originate in such fountain outflows (Fig.[3]and 
Fig. [4}. 

2. We find that the column density profiles of various ions and of 

neutral hydrogen in our simulations are an exponential function of 
the impact parameter (see Section [43] ). Ions with higher ionization 
energy have more extended profiles with the scale height of the 
exponential distributi on tig htly correlated with the ion ionization 
energy (see Fig. H eqs. |8|9| ): h s oc At z « 0, the scale height 

of warm gas traced by low-ionization species, such as Mg II and 
CIV, range in « 0.03 — 0.07itLir, while higher ionization species, 
such as O VI and Ne VIII, have scale heights of « 0.32 — 0.45i? v ir 
(Tableland Fig. [9j. 

3. The predicted correlation of scale height with ionization en¬ 
ergy of ions is in qualitative and quantitative agreement with obser¬ 
vations. In particular, one of our simulations reproduces the radii 
where observations show a sharp turnover in the column density 
distribution for different ions (Fig. [T3] [T4] and |T5j . We emphasize, 
however, that this turnover is not due to a sharp boundary, but is a 
manifestation of the continuous exponential column density profile. 
(Fig. [3}. 

4. We find that the total gas column density profile can be ap¬ 
proximated by a double exponential profile. The outer exponential 
distribution is well traced by the distribution of all ions, although 
column densities of low ions at these radii are predicted to be much 
below current sensitivity limits. However, the model predicts that 
the observed column densities of higher ions track the shape of the 
total gas column density profile at large r quite closely. 

5. We find that the physical extent of the CGM in our simu¬ 
lations evolves slower than the virial radius at z ^ 2. We show 
that column density profiles, in the simulations that match the ob¬ 
served CGM, evolve similarly to the halo scale radius, r s . Thus, 
column density profiles of galaxies at different redshifts can be 
rescaled using r s corresponding to their halo mass. This reveals 
a self-similarity in the simulated and observed CGM across four 
decades of stellar mass and 11 billion years in cosmic time (Fig [7] 

and Fig. ED- 

6. We show that the addition of supernova-produced cosmic ray 
fluid and associated pressure on the gas produces a CGM with pro¬ 
files close to low -z observations. The CGM in this case is much 
more uniform and less patchy compared to the simulations without 
cosmic rays. These differences can potentially be tested by future 
observational measurements of the covering fraction profiles. Over¬ 
all, we find that CR-driven winds in the galactic halo contain cooler 
gas (T < 10 5 K) compared to winds driven by supernova feedback 
without cosmic rays in agreement with previous studies (Fig. 
and[T 6 ]>. 

All but one simulation presented here reproduce the star for¬ 
mation history expected for a typical « L* galaxy reasonably well. 
In agreement with other recent studies, we show that properties of 
the CGM are quite sensitive to the details of the star formation- 
feedback loop. For example, the fiducial simulation of |Agertz“~&| 
|Kravtsov| < |2015[ [20l6] >, ALL_Efb_e010, which produces a very 
realistic central galaxy with the correct stellar mass, size, angular 
momentum, rotation curve, bulge-to-disk ratio, stellar and gas sur¬ 
face density profiles fails to reproduce existing observations of the 
CGM. At the same time, variations of star formation efficiency, en¬ 
ergy per supernova, or introducing cosmic ray feedback appear to 
bring the predicted CGM in better agreement with observations. 
This illustrates that the properties of galaxies and the properties of 
their CGM provide strong complementary constraints on the pro¬ 
cesses governing galaxy formation. Our results clearly show that 


future tests of galaxy formation physics should make use of the 
growing data set of the CGM measurements. Potentially, informa¬ 
tion on the hot gas component via the highest ions or direct X-ray 
imaging could also provide valuable constraints on the models of 
star formation and feedback. 
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